close all
clear all
clc

delta=0.083;
alpha=0.30;
tauK=0.35;
tauI=0.0549;
betta=0.96;
N=0.3;
K=1;
i=1;
omegaN=2.67;
TauN=0; 

%% First Case theta=0.05

theta=0.05;

for tauD=0:0.01:1

treshold=delta/((1-tauD)*(1+tauI));
    
if theta<treshold
    
K(i)=N*((alpha*(1-tauK))/(delta*(((1/betta)-1)/((1-tauD)*theta)+(1+tauI))))^(1/(1-alpha));

K_N(i)=K(i)/N;
p(i)=(delta/(theta))*K(i);
q(i)=delta/(theta);
I(i)=theta*q(i)*K(i);

Taud(i)=tauD;

W(i)=(1-alpha)*(K_N(i))^alpha;
C(i)=(1-TauN)*W(i)/omegaN;
D(i)=(1-tauD)*((1-tauK)*alpha*(K_N(i))^alpha-delta*(1+tauI)*K_N(i))*N;
U(i)=log(C(i))-omegaN*N;
Welfare(i)=U(i)/(1-betta);

else
    
K(i)=N*(((((1+tauI)^(-1))*alpha*(1-tauK))/((1/betta)-1+delta))^(1/(1-alpha)));

K_N(i)=K(i)/N;
p(i)=(1-tauD)*(1+tauI)*K(i);
q(i)=(1-tauD)*(1+tauI);

Taud(i)=tauD;

W(i)=(1-alpha)*(K_N(i))^alpha;
C(i)=(1-TauN)*W(i)/omegaN;
D(i)=(1-tauD)*((1-tauK)*alpha*(K_N(i))^alpha-(1+tauI)*delta*K_N(i))*N;
U(i)=log(C(i))-omegaN*N;
Welfare(i)=U(i)/(1-betta);

end


 i=i+1;

end

figure(1)
subplot(3,5,1)
plot(Taud,K/N)
xlabel('Dividend Tax')
ylabel('K/N')
xlim([0 0.5])

subplot(3,5,2)
plot(Taud,q)
xlabel('Dividend Tax')
ylabel('q')
xlim([0 0.5])

subplot(3,5,3)
plot(Taud,p)
xlabel('Dividend Tax')
ylabel('p')
xlim([0 0.5])

title('\theta=0.05');

subplot(3,5,4)
plot(Taud,D)
xlabel('Dividend Tax')
ylabel('Dividends')
xlim([0 0.5])

subplot(3,5,5)
plot(Taud,Welfare)
xlabel('Dividend Tax')
ylabel('Welfare')
xlim([0 0.5])


%% Second Case - theta=0.11

theta=0.11;
i=1;

for tauD=0:0.01:1

treshold=delta/((1-tauD)*(1+tauI));
    
if theta<treshold
    
K(i)=N*((alpha*(1-tauK))/(delta*(((1/betta)-1)/((1-tauD)*theta)+(1+tauI))))^(1/(1-alpha));
K_N(i)=K(i)/N;
p(i)=(delta/(theta))*K(i);
q(i)=delta/(theta);
I(i)=theta*q(i)*K(i);

Taud(i)=tauD;

W(i)=(1-alpha)*(K_N(i))^alpha;
C(i)=(1-TauN)*W(i)/omegaN;
D(i)=(1-tauD)*((1-tauK)*alpha*(K_N(i))^alpha-delta*(1+tauI)*K_N(i))*N;
U(i)=log(C(i))-omegaN*N;
Welfare(i)=U(i)/(1-betta);

else
    
K(i)=N*(((((1+tauI)^(-1))*alpha*(1-tauK))/((1/betta)-1+delta))^(1/(1-alpha)));
K_N(i)=K(i)/N;
p(i)=(1-tauD)*(1+tauI)*K(i);
q(i)=(1-tauD)*(1+tauI);

Taud(i)=tauD;

W(i)=(1-alpha)*(K_N(i))^alpha;
C(i)=(1-TauN)*W(i)/omegaN;
D(i)=(1-tauD)*((1-tauK)*alpha*(K_N(i))^alpha-(1+tauI)*delta*K_N(i))*N;
U(i)=log(C(i))-omegaN*N;
Welfare(i)=U(i)/(1-betta);

end


 i=i+1;

end

figure(1)
subplot(3,5,6)
plot(Taud,K/N)
xlabel('Dividend Tax')
ylabel('K/N')
xlim([0 0.5])
subplot(3,5,7)
plot(Taud,q)
xlabel('Dividend Tax')
ylabel('q')
xlim([0 0.5])
subplot(3,5,8)
plot(Taud,p)
xlabel('Dividend Tax')
ylabel('p')
xlim([0 0.5])

title('\theta=0.11');

subplot(3,5,9)
plot(Taud,D)
xlabel('Dividend Tax')
ylabel('Dividends')
xlim([0 0.5])
subplot(3,5,10)
plot(Taud,Welfare)
xlabel('Dividend Tax')
ylabel('Welfare')
xlim([0 0.5])


%% Third Case - theta=0.30

theta=0.30;
i=1;

for tauD=0:0.01:1

treshold=delta/((1-tauD)*(1+tauI));
    
if theta<treshold
    

K(i)=N*((alpha*(1-tauK))/(delta*(((1/betta)-1)/((1-tauD)*theta)+(1+tauI))))^(1/(1-alpha));
K_N(i)=K(i)/N;
p(i)=(delta/(theta))*K(i);
q(i)=delta/(theta);
I(i)=theta*q(i)*K(i);

Taud(i)=tauD;

W(i)=(1-alpha)*(K_N(i))^alpha;
C(i)=(1-TauN)*W(i)/omegaN;
D(i)=(1-tauD)*((1-tauK)*alpha*(K_N(i))^alpha-delta*(1+tauI)*K_N(i))*N;

U(i)=log(C(i))-omegaN*N;
Welfare(i)=U(i)/(1-betta);

else
    
K(i)=N*(((((1+tauI)^(-1))*alpha*(1-tauK))/((1/betta)-1+delta))^(1/(1-alpha)));


K_N(i)=K(i)/N;
p(i)=(1-tauD)*(1+tauI)*K(i);
q(i)=(1-tauD)*(1+tauI);

Taud(i)=tauD;

W(i)=(1-alpha)*(K_N(i))^alpha;
C(i)=(1-TauN)*W(i)/omegaN;
D(i)=(1-tauD)*((1-tauK)*alpha*(K_N(i))^alpha-(1+tauI)*delta*K_N(i))*N;

U(i)=log(C(i))-omegaN*N;
Welfare(i)=U(i)/(1-betta);
end


 i=i+1;

end

figure(1)
subplot(3,5,11)
plot(Taud,K/N)
xlabel('Dividend Tax')
ylabel('K/N')
xlim([0 0.5])
subplot(3,5,12)
plot(Taud,q)
xlabel('Dividend Tax')
ylabel('q')
xlim([0 0.5])
subplot(3,5,13)
plot(Taud,p)
xlabel('Dividend Tax')
ylabel('p')
xlim([0 0.5])

title('\theta=0.30');

subplot(3,5,14)
plot(Taud,D)
xlabel('Dividend Tax')
ylabel('Dividends')
xlim([0 0.5])
subplot(3,5,15)
plot(Taud,Welfare)
xlabel('Dividend Tax')
ylabel('Welfare')
xlim([0 0.5])




